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Abstract. We study the fourth order normalized cumulant of height fluctuations 
governed by 1 -|- 1 dimensional Kardar-Parisi-Zhang (KPZ) equation for a growing 
surface. Following a diagrammatic renormalization scheme, we evaluate the kurtosis 
Q from the connected diagrams leading to the value Q = 0.1523 in the large-scale 
long-time limit. 
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1. Introduction 

Growth of a surface or interface has been one of the most important and well-studied 
helds in nonequilibrium statistical physics since a long time Kardar, Parisi 

and Zhang [6] hrst proposed a paradigmatic nonlinear equation for local surface growth 
capable of describing many growth phenomena. The equation, called the KPZ equation, 
is expressed as 



where /i(x, t) is the held of height huctuations, uq is the surface tension that relaxes 


particles from local maxima to local minima, and Aq is the strength of local interaction. 
Here ? 7 (x, f) is the deposition noise with zero average, (^/(x, t)) = 0, and its covariance 
is modeled as a short range correlation 


(77(x, t) t')) = 2Do h'^(x - x') 6{t- T), 


( 2 ) 


with d the substrate dimension. 

The roughness of a self-affine surface is characterized by the width w of the interface 
(or standard deviation w of the height huctuations), given by the dynamic scaling 
relation 



( 3 ) 
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as suggested by Family and Vicsek [S], where L is the size of the interface and /(■) 
is a universal function having asymptotics such that w{L, t) ~ t^ when t ^ and 
w{L,t) ~ for t^ . Here the exponent y characterizes the roughness of the surface, 
z is the dynamic exponent, and the ratio /3 = ^ is known as the growth exponent. The 
roughness exponent y is an important parameter in experiments; adsorption, catalysis 
|7] and optical properties [8] of a thin him are affected by the roughness of the surface. 
The exponents are related via the scaling relation y + ^ = 2 p [ini [m, which is 
independent of the substrate dimension. 

There are various growth phenomena that are believed to be in the KPZ universality 
class on the basis of the numerical values of the scaling exponents pQ. A few of 
them are thin him deposition |T2], bacteria colony growth Haim, huid how in porous 
media |I5], turbulent liquid crystal daE], one dimensional polynuclear growth (PNG) 
[m ca Ea E], slow combustion of a sheet of paper |22[ I23j . In addition, many 
problems are equivalent to the KPZ equation, e.g.. Burgers equation [2l| describes 
the vorticity free velocity, directed polymer in random media [251 126] and in random 
potentials (DPRP) [27], sequence alignment of gene or protein [281129], heat equation 
of multiplicative noise obtained via the Cole-Hopf transformation of the KPZ equation. 

A renormalization group (RG) analysis [6] of the 1 + 1 dimensional KPZ equation 
yielded roughness exponent X = | and dynamic exponent 2 : = | which are 

consistent with various numerical models, e.g., ballistic deposition PP, Eden model 
[301 [311 ESI [33], restricted solid on solid model (RSOS) [1], single step model (SSM) 
P ED- RG calculation in 2 + 1 dimensions is unable to yield the exponents of 
KPZ equation. A number of analytical techniques have been employed to study the 
higher dimensional scaling exponents of KPZ equation, e.g., the mode coupling [35l [36] , 
the operator product expansion [37], the self-consistent expansion [38] and the non- 
perturbative RG [391 HO] . 

Although many physical problems are described by the 1 + 1 dimensional KPZ 
equation, enough attention has not been payed to understand the statistical probability 
distribution of the height fluctuations as it is a challenging task to obtain it directly 
from the KPZ equation. Nevertheless, it is a highly desirable objective and it demands 
extensive theoretical and experimental studies [4]. Prahofer and Spohn [21] studied the 
PNG model for three different initial conditions namely flat, droplet and stationary self 
similar which are in the KPZ universality class on the basis of scaling exponents. For the 
droplet initial condition the distribution was found to be the Gaussian unitary ensemble 
(GUE) Tracy-Widom (TW) whereas for the flat initial condition it is the Gaussian 
orthogonal ensemble (GOE) Tracy-Widom (TW) distribution. The height distribution 
was modeled through the relation h{x, t) ~ Voot+{Tty^^(l) (with the parameter F = 
where 0 is a random variable determined by appropriate random matrices and v^o is the 
rate of growth at long times [T6] . 

Imamura and Sasamoto [H] considered a Brownian motion as an initial condition 
from both sides of the substrate and used the Bethe ansatz and a replica trick to hnd the 
height distribution in terms of a Fredholm determinant. Takeuchi [32] obtained scaling 
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functions smoothly connecting the crossover between the GOE-TW and Baik-Rains Fq 
distributions. He observed that the moments pass through minimum values, known as 
the Takeuchi minima. Such behavior was also noticed via TLC experimental studies 
apart from the numerical study of the PNG model. 

In spite of the same scaling exponents, the statistical behavior of stochastic 
processes may be different due to difference in the underlying probability distribution 
functions (pdf) that incorporate the basic features of a dynamical process [13]. In 
principle, the pdf can be calculated analytically by solving the Fokker-Planck equation 
[mils] corresponding to the KPZ equation. However, this is practically infeasible due to 
the nonlinear term. As an alternative, a few higher order moments can be calculated to 
understand the statistical behavior of the interface and the corresponding universality 
class. It may be noted that the measurement accuracy of higher order moments is 
greater than that of the scaling exponents |16] in experiments. 

In this paper, we consider the 1 + 1 dimensional KPZ equation with a flat substrate 
and focus on the evaluation of the fourth order cumulant relevant to the kurtosis in 
stationary state. We apply a diagrammatic approach to find a renormalized expression 
for the loop integral corresponding to the fourth order cumulant in the large scale long 
time limit. 

The paper is organized as follows. Section 2 defines the moments and cumulants 
and states the relations between them. Section 3 is devoted to the calculation of fourth 
order cumulant. In Section 4, the calculation of excess kurtosis is presented. Section 
5 presents a discussion and conclusion and a comparison the kurtosis value with other 
numerical and experimental findings. 

2. Moments and Cumulants 

Higher order moments and cumulants are usually defined in terms of a generating 
function. The moment generating function Z{fl) is defined as 

°° on 

Z(/?) = (e'’'>) = ^+(ft">, (4) 

n\ 

n=0 

where (/i") is the nth order moment of a random variable h, given by 

/ +00 

h^P{h) dh, (5) 

-OO 

for a given normalized probability distribution function P{h). The cumulant generating 
function F{(3) is defined as 

°° on 

F(/j)= In Z(/J) = 5;+{/!”)„ (6) 

n=l 

where {h"')c is the nth order cumulant. The moments are related to the cumulants as 
(h) =(h), 

{h^) = {h^)c+{h)l 
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{h^) — {h?)c + ‘i{h)c{h^)c + (^)c 

(fi^) = (h^), + 4{hUh% + 3{hYc + Q{h)l{h\ + {h)t. (7) 

The fluctuating interface field /i(x, f) has a zero mean; (/i(x, f)) = (/i(x, f))c = 0. Hence 
the second and fourth order cumulants [271143] are related to the moments as 

{h^)c={h^), (8) 

and 


{h^).= {h^)-3{hY. 

Kurtosis Q is defined by the “normalized” fourth order cumulant as 

{h% {h‘) 


0 = 


- 3. 


(9) 


( 10 ) 


(h?)i {h^y 

In the diagrammatic approach, a connected diagram with n external legs corresponds 
to the nth order cumulant. Henceforth, we focus on the fourth cumulant instead of the 
moment and therefore we evaluate the connected loop diagram to obtain the cumulant 
in order to evaluate the kurtosis. 


3. The Fourth-order Cumulant 

The fourth order cumulant (/i^(x, t))c of the height fluctuations measures the degree of 
flatness of the probability distribution. Using the Fourier transform 

the KPZ equation ([T]) is written in momentum and frequency space as 

{-iu+vofl^) h%u) = 77(k,a;)-y j [q-(k-q)] /i(q, H) /i(k-q,cn-H). 
This equation forms the basis of developing a perturbation theory. 



Figure 1. Feynman Diagram corresponding to the fourth cumulant 


The expression for the fourth cumulant is given by the connected fourth order 
correlation in momentum and frequency space as 


{/t^(x,t))c = 


f f d^+^k f 


(27r; 


d+1 


(2y 


d+1 


(2^: 


d+l 


(2^: 


d+l 




( 13 ) 
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where x = (x, t) and ki = (kjjCJj) are vectors in (d + l)-vector notation. We use the 
diagrammatic approach and obtain a connected loop diagram for the fourth cumulant 
as shown in Fig. 1, where a wiggly line represents the correlation and a solid line the 
response function. We express the correlation in accordance with the connected loop 
diagram in Fig. 1, so that 






J { 2 nY+^ J { 27 vY+^ J ( 27 r 
Lflki, ^ 2 , kf) G{—ki — k2 — kflj, 


gd+il ^ ^ ^ 

- ^G{k,)G{k2)Giks) 


(14) 


where L 4 represents the renormalized loop corresponding to the amputated part of the 
diagram (excluding the external legs). Its bare value is written as 

L®(fci,^2,fc3) = 16 (2Do)^y (q-ki)][q-(k 2 + q)] 

[(q + k 2 ) • (ks + k 2 + q)] [(k4 + ka + k 2 + q) ■ (q + k 2 + ks)] 
GoiQ)Go{~Q + ki)GQ{—q)Go{q + ^2)G'o(~? ~ k2)Go{q + k^ + ^ 2 ) 
Go{q + ^2 + ^3 + k4)Go(—q — k2 — k^), ( 15 ) 


where the prefactor 16 is a combinatorial factor. 

Now we evaluate Eq. [T5]by performing frequency and momentum integrations where 
the momentum is restricted to the shell Age"'’ < q < Aq. This leads to 

y 


c; r)4 

T<(r']- 0 0 


( 16 ) 


where Ka = with Sd the surface area of a unit sphere in d dimensional space. 

We follow Yakhot and Orszag’s [17108] scheme of renormalization without rescaling 
and obtain from Eq. [12] the differential equation 

dL^ 5 X^oD\r) 


dr 2 tt u'^{r)A^{r)' 
in one dimension. 

Noting that A(r) = Aoe“'’ and using the asymptotic expressions 


(17) 


z/(r) = Ac 


Dn 


27ri/oAo 


y /2 


and 


D{r) = 


XqDq 


Dn 


y /2 


z/Q V 27rz/oAo 
obtained from renormalization-group calculations 
differential equation flTTl) as 


Lflr) = 


IOtt / 2D^ 


\7rnlAlJ 


Jr /2 


(18) 


(19) 


we obtain the solution to the 


( 20 ) 
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The factor is interpreted as a momentum ki in the large scale limit. Since the 

loop integral must be symmetric with respect to interchange of the external momenta 
ki, k 2 and /cs, we replace Aoe“^ by the fully symmetric form k\^^ky^ky^. Thus we have 


L4(ki, 0; k2, 0; ka, 0) — 


IOtt 

~Y 


Ar 


9 r)5 \ 

\ ,^- 7 / 6 ,- 7 / 6 ,- 7/6 


VTl^n 


Ail ^2 ^3 


( 21 ) 


when the external frequencies are zero. To obtain the frequency dependency, we take 
the scaling function as 


bJi 


Ph. 1 p I = 


( 22 ) 


where the renormalized response function is given by = [—iui + v{ki)\ffl\~^ 

with the renormalized surface tension 


I'iki) = Ao 




27ri/n 


(23) 


The above scaling relation fl2^ obeys consistency with the zero frequency limit and with 
the real valuedness of Age”’', in addition to being of the correct dimension. Now with 
the aid of Eq. ([22]), Eq. fl2l]l is modihed to the frequency dependent form 


Lflh,k2,k) = ^Xo 


9 7~)^ \ ^ A 

,+ 7 / 6 , 17 / 6 , 17 / 6 , 2 , 


TTZ/n 




„"(fc,)|G(l2,a,2)|"„"(fc3)|G(l3.ws)^ 


(24) 


representing the renormalized loop diagram in Fig. 1. We substitute the renormalized 
expression for the loop diagram from Eq. [21] in Eq. HI and treat the external legs 
representing the response functions as renormalized. Carring out the frequency 
integration, we obtain 


where 


{/i^(x,t))c 


10 

y 



r*oo poo 

dko / 


dk^^{ki,k2,kfl 


^{ki,k2,kfl 


ujki, k 2 , ks) 

v{ki,k2,k3) 


(25) 


(26) 


with 


u{ki, k 2 , ks) = (3(|kyA + + [kgl^A) 5(|kyA + |k 3 |^A)|k^ + ks + 

+ 7(|ky + |k2p + IksH |ki + k2 + kgp/^ + 17(|ky |k2|3/2 ^ |ki|3A|k2|" 
+ Ikyiksi"/' + |ky/2|]^3|" + Ika^kgl"/' + |k2|"/>3|") 

+ 22(|kyA|k2|"/' + |kyA|k3|3A + |k2|3A|k3|3A)|ki + ks + 

+ 90|kyA|k2|3A|k3|3A + |ki + ks + kgl^A) (27) 

v{k,, k 2 , ks) = 64|ki|A6|k2|A6|k3|7/6(|k3|3/2 + |k2|3A + |k3|3A + |ki+k2+k3p/2)A(28) 
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The symmetry of the function ^{ki, ^ 3 ) allows us to write Eq. EBlas 

flOO flOO 

dki dk2 dks [2^{ki, k2, ks) + Q^{-ki,k2,k3)] (29) 
J ^ J fL 

where jj, is an infrared cut off. Now we write the integrations separately as 


{h (x,f))c = — 


10 / D, 


7 V 2ni>(] 


A(+) 



dk2 dk3^{ki,k2,k3) 


and 




dki / dk2 dfcs <h(—A:i,/c2, ^3)- 


These integrals are expected to be of the forms 

= ai/r^ 

^2(/() — 


(30) 

(31) 

(32) 

(33) 


where ai and 02 are dimensionless constants. Substituting Eqs. |32] and |33] in Eq. [291 
yields 

10 / D \^ 1 

(+x,t)),= -[2a,+6a2l(^) ^ (34) 

We evaluate the integrals in Eqs. |30]and[3T]by integrating over ki, k 2 and k 3 numerically 
and obtain the dimensionless constants as 


and 


ai = lim [/i^Ji(p)] = 0.007505 


(35) 


02 = lim |y/ 2 (A)l = 0.026297, (36) 

/4^0+ 

the numerical values having converged to the above values for decreasing values of the 
parameter /i very close to zero. 


4. The Kurtosis 


Having calculated the fourth-order cumulant given by Eq. fl34p . we now need the value 
of the second-order cumulant to determine the value of kurtosis given by Eq. flTOl) . The 
second moment is written as 



d‘^k2 du2 

{2'K)d+l 


(h(ki, a;i)h(k2, wa)) 


Homogeneity in space and time implies that the correlation takes the form 


(/i(ki,(Ui) h(k2,(Ua)) — Q(ki, uji) (27r)'^5^(ki + ka) (27r)h(a;i + cua)- (38) 


The correlation function Q(k, uj) can be written in accordance with the one-loop diagram 
in Fig. 2, so that 


(/i^(x,t))c 



d‘^k du 

(27rY+i 


G(k,u)L2(k,u)G{-k,-u) 


(39) 
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Figure 2. Feynman diagram corresponding to the second moment 


The renormalized value of L 2 , given by the loop integral coming from the amputated 
part of Fig. 2, can be obtained as 

L2(k,c) = ^4|G(k,i..)|, (40) 

by means of a renormalization scheme as shown in Ref. [19]. This yields the second 
cumulant as 


{/t^(x,t))e = - 


TT 


Do 

27rzyn 


1 

h 


Employing Eqs. (I25ll and (HTjl in the dehnition of kurtosis leads to 


Q = 




(oi + 802) 


(41) 


(42) 


(/l2(x,t))2 

free from the infra-red cutoff /i and the model parameters uq, Dq, and Aq. Employing 
the numerical values of ai and 02 from Eqs. [35] and [36| yields the kurtosis as 


Q = 0.152267. 


(43) 


We compare this value with other stationary values of kurtosis in Table 1. 


Table 1. Values of Kurtosis in (1 + 1) dimensions. 


System of study 

Method 

Kurtosis 

Reference 

PNG (stationary) 

Numerical 

0.289 

1211 

gl DPRM (stationary) 

Numerical 

0.278 

ra 

KPZ (present calculation) 

Analytical 

0.1523 

Eq. (|13D 


5. Discussion and Conclusion 

In this work, we followed a perturbative renormalization approach to evaluate the 
Feynman diagram in Fig. 1 that represents the fourth cumulant, noting that contribution 
to cumulants are given by the connected loop diagrams. To calculate the connected 
diagram, we hrst renormalized its amputated part in the limit of internal momenta 
and frequencies much greater than the external ones and followed Yakhot and Orszag’s 
iterative renormalization scheme without rescaling |18]. This results in a scale dependent 
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function for the amputated part of the loop diagram at zero frequency. Introducing a 
frequency dependent scaling function that smoothly joins with the obtained scaling 
behavior and that preserves the property of real valuedness, we evaluated the resulting 
frequency and momentum integrations to obtain an expression for the fourth order 
cumulant depending on the infrared cutoff in the momentum integration. The evaluation 
of the second cumulant is simpler which also depends on the the infrared cutoff as shown 
in Ref. |19]. The behaviors of both the fourth and second order cumulants are seen to 
be exactly in accordance with the well-known behavior of the nth moment, namely 
Wn ~ where L is the substrate size (so that /r ~ L~^) and y is the roughness 
exponent. These cumulants yield the value of kurtosis Q = 0.1523 which is free from 
the infrared cutoff p and model dependent parameters Dq , i/q, and Aq, due to their 
exact cancellations. 

Incidentally, it may be worth mentioning that the calculation of cumulant 
amplitudes have been regarded as an important objective by some researchers. For 
example, Krug et ah 1271 studied the growth of an interface via simulation of the 
single step model with flat initial condition and obtained the estimates for the cumulant 
amplitudes as C 2 = 0.404 ± 0.013 and C 4 = 0.020 ± 0.002, suggesting a kurtosis value 
Q = % = 0.123 ± 0.020. Furthermore, they showed that C 2 is higher for the steady 
state case, namely, C 2 = 0.712 ± 0.003, due to the additional role of fluctuations in 
the initial conditions besides those operating during the growth. Earlier, Hwa and 
Frey [51] had obtained C 2 = 0.69 via mode-coupling calculations. Tang [52] estimated 
C 2 = 0.725 ± 0.005 using a Monte-Carlo simulation with the single-step model, giving 
an excellent estimate for the Baik-Rains constant [SUES]. A slightly different result, 
C 2 = 0.71, was obtained through mode-coupling calculation by Amar and Family [54] . 

In Table 1, our calculated kurtosis value is compared with other stationary values. 
Prahofer and Spohn [21] studied the PNG model belonging to the 1-1-1 dimensional 
KPZ universality class to hnd the effect of initial conditions on the statistical behavior 
of growing interfaces. They obtained different kurtosis values for different initial 
conditions, namely, Q = 0.0934 for curved, Q = 0.1652 for flat, and Q = 0.289 for 
stationary initial conditions. Halpin-Healy and Lin [50] studied the DPRM with different 
conhgurations such as point-to-point, point-to-line, and stationary cases, leading to 
GUE-TW, GOE-TW and Baik-Rains Fq distributions, respectively. Through gi DPRM, 
they estimated the stationary kurtosis value to be Q = 0.278. 

Apart from the values shown in Table 1, there are a few studies where the full 
Baik-Rains distribution is obtained in the context of KPZ stationary state. Takeuchi 
[42] obtained, via a numerical simulation of the PNG model and an experiment on 
the TLG, a universal function that undergoes a crossover from a transient state to the 
stationary regime. Miettinen et al. [23] carried out an experimental study of the slow 
propagation of the combustion front on a sheet of paper to investigate upon the front 
fluctuation distribution. Their experimental data for the transient and stationary states 
were found to ht well with the GOE-TW and Fq distributions, respectively. On the other 
hand, Halpin-Healy and Lin [50], studied the distributions of deposition models (BD, 
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SSM, and RSOS) that strongly agree with DPRM/SHE results and they turn out to 
be the Baik-Rains distribution. These investigations establish that the kurtosis value 
in the stationary state must be the same as the universal kurtosis of the Baik-Rains 
distribution (for example, Q = 0.289 [2Tj or 0.278 [50]). 

In our calculation based on the RG scale elimination scheme, we obtained the 
renormalized quantities in the large-scale and long-time limits. Thus our methodology 
inexorably selects the stationary regime. However, the resulting kurtosis value does not 
agree well with the stationary kurtosis value and it is distinctly lower than the Baik- 
Rains value, Q = 0.28916. In our simplihed scheme of calculation, we obtained the 
fourth cumulant at one-loop order. While a one-loop scheme for the third cumulant 
was nearly successful in estimating the stationary skewness value |19|, a higher order 
calculation would appear to be more appropriate for the fourth cumulant. Noting that 
the perturbation expansion is about a Gaussian state, the Gaussian behavior seems to 
play some role with the increase in order of the cumulant, lowering the estimated value. 
This opens the door for further investigations into the rather unexplored stationary 
KPZ problem that is expected to illuminate upon its semantic relation with Baik-Rains 
distribution. 

We conclude by recalling that there are many growth processes governed by 
nonequilibrium dynamics that are believed to be in the KPZ universality class on 
the basis of scaling exponents. However, as shown by the PNG and TLG studies, 
the statistical behavior, that is the probability distribution, depends on the initial 
conditions. The full probability distribution function of the KPZ height fluctuations has 
never been studied analytically due to the inherent difficulty in the problem. However, 
relevant information about the probability distribution can be obtained through the 
study a few higher order moments and cumulants. Thus the present theoretical study 
may be viewed as an initial attempt at the quantification of higher order statistical 
property of the 1-1-1 dimensional KPZ dynamics. We hope that such analytical scheme 
would be useful for the study of statistical behavior of other important stochastic 
processes as well. 
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